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Preface 


This report is a part of the research and development project on air- 
craft de-icing by the electromac .etic impulse method. This project has been 
sponsored by the Lewis Research Center of the National Aeronautics and Space 
Administration under grant number NAG-3-284. The grant administrators were 
Mr. John J. Reinmann and Dr. R. Joseph Shaw. For the previous five years, 
many tests had been run in the NASA Icing Research Tunnel and three sets of 
flight tests were performed. These, plus various laboratory tests, have 
resulted in a semi-empirical technology for designing an electro- impulse de- 
icing (EIDI) system. 

However, the empirical method is inadequate when a very different 
geometry, material or size is encountered. A computative solution is needed 
which permits prediction of the de-icing effect for a given configuration 
and electrical circuitry. This report, which is principally the Ph.D. 
dissertation of Robert A Henderson under the direction of Prof. Robert L. 
Schrag, attempts to do the first part of a full computer simulation of EIDI. 
The pressure/time produced by the method in the report would be necessary 
input for a computer code giving the structural dynamic response of a given 
configuration. The configurations in mind are leading edge portions of 
aircraft wings, engine nacelle inlets and rotor blades. Applications, 
however, are not limited to these uses. 

The authors acknowledge the assistance of NASA-Lewis Research Center, 
both for the support of the whole EIDI project at Wichita State University 
and for the opportunity for Dr. Henderson to work at NASA-Lewis during the 
summers of 1985 and 1986. The assistance of Drs. R. Joseph Shaw, Bill Ford 
and Avram Sidi during that time is gratefully expressed. 
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CHAPTER ONE 


INTRODUCTION 

1.1 ET.FTrrRO-TMPULSE PHENOMENON AND DE-ICING 

In a paper presented January 8, 1986, at the 24th AIAA Aerospace 
Sciences Meeting in Reno, Nevada, [1] author R. D. Rudich cited aircraft 
icing as the direct cause of four of the thirty two weather related fatal 
accidents reported in the paper. While this may not seem like a signifi- 
cant number, the loss of human lives in these four accidents could have 
been prevented if the aircraft involved had been able to cope with the 
icing conditions they encountered. 

The purpose of this report is to present methods of analyzing the 
electromagnetic aspects of a new method of de-icing both private and 
commercial aircraft. This new de-icing method, referred to as Electro 
Impulse De-Icing (abbreviated EIDI), holds the promise of being superior 
to present aircraft de-icing methods in terms of the energy expended in 
removing accreted ice [42]. 

The EIDI concept is not new. In May 1939, Great Britain issued a 
patent to Mr. Rudolf Goldschmidt covering the basic EIDI mechanism [2]. 
However, no commercial development of an EIDI system in the free world 
proceeded from this patent. It has only been within the last 5 years 
that the unavailability of bleed air from the high bypass ratio engines 
used on the next generation commercial aircraft has caused attention to 
be directed to the use of an EIDI system for removing ice. 

The simplest practical EIDI system consists of a spirally wound coil 
of rectangular cross section conductor mounted with its axis of symmetry 
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perpendicular to the metal surface to be de-iced. An initially charged 
capacitor is discharged through this coil, and the resulting magnetic 
field from the coil's current causes eddy currents in the metal. The 
force exerted on these eddy currents by the coil's magnetic field is 
initially in such a direction as to cause the coil and the metal surface 
to separate. It is this force that causes ice on the metal surface to 
crack, and subsequently to be removed from the surface. 

A considerable body of literature concerned with the electromagnetic 
aspects of a coil placed next to a conducting surface has accumulated. 
Levy [8] and Grover [9] present methods of calculating terminal imped- 
ances. Dodd and Deeds [10] discuss both impedance calculations and eddy 
current distributions. Onoe [11] is apparently the first researcher to 
apply a Hankel transformation in the calculation of impedances. Onoe's 
method is further developed and extended by El-Markabi and Freeman [4] to 
include calculation of the force between the coil and conductor when the 
coil current is a sinusoid. Bowley et. al. [43] discuss the use of the 
magnetic vector potential in calculating the impulse delivered to the 
conductor by a single cycle of exponentially damped sine wave current in 
the coil. Lewis [44] summarizes the use of the Bowley method for design- 
ing a coil to deliver a specific inpulse in an electro- inpulse de-icing 
installation. 

1.2 SCOPE OF REPORT 

An experimental set-up of a prototype EIDI system was assembled at 
The Wichita State University by Dr. Robert Schrag. An experimental study 
of the electro-impulse phenomenon in this system was made by Dr. Schrag 
utilizing field diagnostics methods [3]. Axial and radial components of 
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the magnetic field were measured on both sides of the rigidly held alumi- 
num "target” plate. The data were then used to deduce the total mechani- 
cal force versus time, and the mechanical impulse strength. Impulse 
strength was also measured directly with a ballistic pendulum. This 
experimental study provided results which will be used to verify theoret- 
ical predictions from the mathematical model used in this report. 

In this report, the physical phenomena involved in the prototype 
EIDI system of Dr. Schrag's experiment will be investigated analytically 
and numerically. Specifically, the following tasks will be undertaken. 

1. A mathematical model for the total electrical problem will be 
devised. It will employ a transmission line analogy to handle the elec- 
tromagnetic field portion of the system, and a frequency domain model for 
the circuit portion. 

2. The math model will be solved by computer for the specific set 
of conditions that existed in Dr. Schrag's diagnostics experiment, and 
the calculated and experimental results will be compared. 


1.3 ORGANIZATION OF REPORT 

Figure 1 summarizes the analysis structure developed in this report 
to predict the behavior of the prototype EIDI system briefly described in 
Section 1.1. Each of the blocks in this figure represents a stage in the 
procedure for determining the force-time profile on a rigid coil placed 
next to a fixed conducting plate (the "target") when a capacitor is dis- 
charged through the coil. 

In Block 1 of Figure 1, Maxwell's equations are written for the coil 
and metal target of the physical system. A simple model of the coil that 
stresses the geometrical symmetry of the system is proposed , eliminating 
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FIGURE 1 

Analysis Flow Diagram 
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several terms from the Maxwell equations. Chapter 3 discusses this 
modeling of the coil and metal target, and shows that application of 
Fourier and Hankel transforms to the model equations results in a concep- 
tual replacement of the field problem with an infinite set of transmis- 
sion line problems (Block 2). 

In Block 3 of Figure 1, the analysis of the transmission line model 
is performed. The methods of this analysis, for steady state sinusoidal 
conditions, may be found in Chapter 4. Chapters Z and 4 provide the 
basic theoretical developments, with the results of the analyses in 
Chapter 4 applied in Chapter 6 to Dr. Schrag's prototype experimental 
EIDI system, described in Chapter 5. Modeling of the circuit used to 
provide energy to the coil is discussed in Chapter 2 for the specific 
experimental system of Chapter 5. 

Once the transmission line model of the coil and target has been 
established in Block 3, Figure 1 shows two possible next steps in the 
analysis. One of these next steps, calculation of the true field values 
(Block 7), requires a knowledge of the coil current. Accordingly, one 
proceeds from the transmission line analysis in Block 3 to the Block 4 
calculation of the coil impedance, required for calculating the coil 
current. Use of the transmission line model for calculating the part of 
the impedance of the coil that is due to the interaction between the coil 
and target is discussed theoretically in Section 4.2 of Chapter 4, and 
applied to the impedance calculation of the coil in the prototype EIDI 
system in Section 6.2 of Chapter 6. 

Coil impedance obtained in Block 4 is added to the calculated skin 
effect loss resistance in the coil. With this total coil impedance, the 
problem of calculating coil current, using the circuit models presented 
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in Chapter 2, is addressed. This is Block 5, discussed in Section 6.3 of 
Chapter 6 for the prototype EIDI system described in Chapter 5. 

Knowledge of the frequency spectrum of the coil current allows the 
calculation of the true field values in Hankel-Fourier space, shown in 
Block 7 and developed theoretically in Section 4.3 of Chapter 4. Per- 
forming an inverse Hankel transformation on these fields results in 
Fourier space fields. An inverse discrete Fourier transform applied to 
the Fourier space fields provides the time domain behavior of the elec- 
tric and magnetic fields, shown in Block 8. This procedure is discussed 
in Section 6.4 of Chapter 6 for the prototype EIDI system. 

Finally, the real space axial and radial magnetic induction fields 
on the coil side face of the metal target provide the necessary informa- 
tion for calculating the total force on the target as a function of time, 
as shown in Block 9. Development of the force equation is in Section 4.5 
of Chapter 4, with Section 6.5 of Chapter 6 describing the numerical 
implementation of this force calculation for the prototype EIDI system. 
Calculation of the inpulse delivered to the target is discussed theoreti- 
cally in Section 4.6 of Chapter 4, and its application to the prototype 
EIDI system is described in Section 6.6 of Chapter 6. 

Chapter 7 contains a summary of the results obtained, and some sug- 
gestions for inproving the procedures for modeling the electrical aspects 
of an EIDI system described in this report. 

A magnetic tape listing of the FORTRAN algorithms developed during 

the research described in this report are available at a nominal charge 

from the authors, at the address below. 

Department of Electrical Engineering 
The Wichita State University, Box 44 
Wichita, Kansas 67208 
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CHAPTER TWO 


THE ELECTRIC CIRCUIT MODEL 


Modeling of the circuit providing power to the de-icing coil in an 
EIDI system is performed with the desired output from the model in mind. 
Since it is the coil current in conjunction with the physical configura- 
tion of the coil and metal target that determiner; the magnetic fields re- 
sponsible for the forces on the target, the current in the circuit was 
chosen as the primary variable. 

Testing of a prototype EIDI system began at The Wichita State Uni- 
versity in 1982. The part of this system that provides power to the coil 
is shown in Figure 5 (page 34). This system was chosen as the basic 
physical configuration for which a circuit model would be derived. The 
circuit model is then analyzed to determine the coil current. A simple 
circuit model that is suggested by the physical system in Figure 5 is 
shown in Figure 17 (page 46). This is the basic circuit model, valid 
when the clamp diode across the capacitor is not conducting. 

Justification for modeling the EIDI system as a lumped parameter 
circuit was provided by the experimental verification of the absence in 
the signals present of any significant frequency components having wave- 
lengths comparable to the physical dimensions of the system. Because of 
the complex electromagnetic interaction between the coil and the target, 
it was felt that a time domain model of the coil's terminal v-i charac- 
teristics would be too complex to be of much use. Consequently, a fre- 
quency domain approach was chosen for analyzing the circuit, making the 
model of the coil a linear frequency dependent impedance. 
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The presence of the SCR and the clamp diode across the capacitor 
makes the circuit nonlinear, so that straightforward Fourier (frequency 
domain) techniques are not applicable. This difficulty is circumvented 
by performing a piecewise linear £ lalysis of the circuit. In this analy- 
sis, the physical circuit is modeled by one of two possible circuits, 
depending on the state of the clamp diode. Initially, when the SCR has 
just been triggered, the diode is assumed off (an open circuit) and the 
model of figure 17 (page 46) is used for analysis. P^th the current and 
the capacitor voltage in this circuit are calculated. When the capacitor 
voltage first becomes negative, the claitp diode paralleling the capacitor 
is modeled as caning into immediate forward conduction. The diode then 
acts as a short circuit, resulting in the circuit model shown in Figure 
18 (page 49). All subsequent (in time) circuit calculations are per- 
formed using this model. 

Theoretical justification for modeling both the SCR and the clamp 
diode as simple on-off switches is now provided. However, the ultimate 
justification for such an outright dismissal of the effects of both the 
SCR and the clamp diode in determining the coil current comes from ob- 
serving how closely the predicted current (using the models that ignore 
the non-ideal nature of the SCR and diode) agrees with the measured 
current. This will be shown in Chapter 6. 

In a practical EIDI system, the voltage on the energy storage capac- 
itor prior to its discharge is 1000 to 1500 volts. When the SCR is trig- 
gered into conduction, its voltage drop is on the order of 1 volt, which 
is small compared to the capacitor voltage, and so may be ignored in 
determining the circuit current. In addition to ignoring the voltage 
drop across the SCR, the dynamics of the SCR are also not modeled. Such 
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dynamics are, in the EIDI circuit, primarily manifested in the failure of 
the SCR to trigger into instantaneous full forward conduction upon appli- 
cation of a forward gate current. However, modem SCR design techniques 
[5], [6], [7] have resulted in tum-on times that are short compared to 
the time required for a significant change in the coil current in the 
experimental EIDI prototype system. The SCR used in the prototype system 
had a di/dt rating of 800 amps/microsecond , while the maximum observed 
rate of change of current in the circuit was 15 amps/microsecond. 

A similar consideration of the voltage levels in the circuit results 
in the conclusion that the approximately 1 volt forward drop across the 
clamp diode is not a primary factor in determining circuit current. When 
the diode is off, its transition capacitance is insignificant compared to 
the 600 microfarad capacitance of the energy storage capacitor in paral- 
lel with it. Reverse leakage current in the diode is ignored due to this 
large energy storage capacitor and the relatively small time in which the 
electrical events of interest take place in the circuit. With the diode 
in forward conduction, the sum of its transition and diffusion capacitan- 
ces are small enough that they also may be ignored in the circuit model. 

In constructing the EIDI prototype experiment configuration, care 
was exercised to minimize parasitics in the circuit. Special low induc- 
tance and resistance cable was used to connect the energy storage capaci- 
tor to the coil. However, both of these cable parameters were measured 
in the prototype system and are taken into account in the circuit model. 
The energy storage capacitor, which was physically several capacitors in 
parallel, used copper strap for wiring connections to minimize inductance 
and resistance parasitics. The equivalent series resistance and equiva- 
lent series inductance of the capacitors were felt to be small, and are 
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not modeled. If, in a particular EIDI installation, these parasitics are 
not small, they may be taken into account by adding lumped elements in 
series with the energy -orage capacitor in the circuit model. 

Accurate frequency domain modeling of the coil is the most diffi- 
cult feat in constructing the circuit model, and is discussed in Chapter 
4 . 
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CHAPTER THREE 


TRANSMISSION LINE MODEL OF THE FIELD PROBLEM 

3.1 INTRODUCTION 

The transmission line model of the coil and metal target is the 
heart of the prototype EIDI system model. It is this model that is ex- 
pected to account for the complex electromagnetic interactions between 
the coil and target. These interactions help establish the impedance 
presented at the coil terminals, and so are a factor in determining the 
coil current. This current is important in determining the force on the 
target. 

3.2 GEOMETRY OF A PROTOTYPE EIDI SYSTEM 

Figure 2 shows the profile of the physical coil-metal plate con- 



figuration to be modeled. The coil has the shape of a short (h«R2) 
thick walled (R 1 «R 2 ) hollow cylinder with an inside radius and an 
outside radius R 2 , whose axis is perpendicular to a flat metal target of 
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thickness d extending to infinity in all radial directions. An air gap 
exists between the coil and target, both of which are assumed rigid and 
stationary in space. 

Most of the coils used in EIDI applications have been round, at 
least prior to a possible bending of the coil to conform to the curved 
leading edge of a wing. The flat geometry of Figure 2 is a reasonable 
model of such a coil, provided that the radius of curvature of the coil, 
after being shaped to conform to the wing, is much greater than the coil 
outer radius R£. 

The assumption that the coil and target remain separated by a fixed 
distance would be reasonable if the initial separation distance was much 
larger than the maximum change in separation distance obtained when the 
capacitor is discharged through the coil. Such an inequality in separa- 
tion distances may not hold in a practical EIDI installation. Alterna- 
tively, the fixed separation distance assumption would be justified if 
the force "impulse" delivered to the plate by the coil was so short that 
the plate acquired only a small velocity, with negligible displacement, 
for the duration of the "inpulse". This does not generally happen in a 
practical EIDI system. In fact, a well designed installation has the 
target move from zero to maximum displacement within the duration of the 
"impulse". With ice loading present, however, the displacement may be 
small corrpared to the initial separation distance. 

3.3 THE FTFTD EQUATIONS 

The approach chosen for modeling the coil and target in the proto- 
type EIDI system is presented in a paper by El-Markabi and Freeman [4]. 
This approach is now described, emphasizing those aspects of the theory 
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that are most appropriate to the analysis of the model of the prototype 
system of Figure 2. The reader is referred to El-Markabi and Freeman 
[4] for the more general theory. 

Current in the coil of Figure 2 is assumed to be entirely phi di- 
rected. A slight modeling error is introduced by this assumption, since 
a radial component of current actually exists in a practical coil due to 
the spiral winding of the ribbon conductor. By neglecting this small 
radial current, a model having azimuthal symmetry is obtained. Conse- 
quently, the model shows no phi dependence in any of the field quanti- 
ties, the electric field contains only a phi component, and the magnetic 
field contains only axial and radial components. 

Because most of the spectral energy of the current in a practical 
EIDI coil is confined to relatively low frequencies, a quasi-static 
situation is assumed. The displacement current term in Ampere's Law is 
ignored. 

With these assumptions. Maxwell's equations written for Figure 2 
become 


9E$(i,r-0 _ 3H r (z, r ,'t^ 

oi e ” r at 

9a 9r 

Partial differentiation with respect to time can be eliminated from these 
equations by taking the Fourier transform. Differentiation with respect 
to time is then replaced with multiplication by the £a> operator. Per- 
forming this transformation leads to the following equations (note that 
new symbols have not been introduced for the transformed field quantities; 
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instead, their time variable argument t has simply been replaced with the 
Fourier transform variable omega) 
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Equation (1) contains only the phi component of E and the radial 
component of H. By combining equations (2) and (3) in such a fasti ion as 
to eliminate the axial component of H, another equation containing only 
the phi component of E and the radial component of H results. This 
elimination yields 

- <r~E. ^ 

= 
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a mathematical tool that is of considerable import- 
This tool is the Hankel transform of order n, de- 


Although the use of an additional transform introduces difficulties 
of its own, it provides an even greater amount of simplification, and 
makes possible the transmission line model of the coil-target electromag- 
netic field problem. 

If a Hankel transform of order 1 is applied with respect to the 
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variable r in equation (4), the following result is obtained [11], [12]. 
Note that new symbols have again not been introduced for the transformed 
quantities; instead, the spatial variable argument r has been replaced 
with the Hankel transform variable lambda. 
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The advantage of equation (5) over equation (4) is that equation (5) 
contains partial derivatives with respect to only the z coordinate, in 
effect becoming an ordinary differential equation (if one allows "con- 
stants of integration" for such an equation to be arbitrary functions of 
all variables except z). If the same Hankel transform is applied to 
equation (1), 

= ^uU r (*,X,co) (*) 

Equations (5) and (6) are two coupled ordinary differential equations 
that are recognizable as the canonical transmission line equations. The 
solutions of these equations are well known [13], [14], [15]. 


3.4 DEVELOPMENT OF THE TRANSMISSION LINE MODEL 

Because equations (5) and (6) are identical in form to the equations 
describing voltage and current on an ordinary transmission line, symbols 
V and I (each of which is a function of position z along the coil axis, 
the Fourier variable omega, and the Hankel variable lambda) are now 
introduced to represent the phi component of E and the radial component 
of H respectively. Then equations (5) and (6) become 
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Differentiating equation (7) with respect to z. 
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Substitution of equation (8) into this last equation yields 
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Similarly, differentiate equation (8) to obtain 
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Substitution of equation (7) into this equation yields 
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Now define the complex propagation constant gamma in Fourier-Hankel space 
as 


y =• -/x a + 

so that the general solutions to equations (9) and (10), assuming complex 
sinusoidal time variation of V and I, may be written as 
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where I 0+ , Iq- and V 0 +, V 0 _ are complex phasor functions of the Hankel 
variable lambda (the carp lex sinusoidal time variation has been suppres- 
sed in these solutions). 

Note that there are an uncountably infinite set of transmission 
lines in the Hankel space model. Each transmission line is associated 
with a different value of lambda, with lambda real and non-negative. The 
physical origin of this infinite number of transmission lines lies in 
having "compressed" the infinite radial variation in the real space field 
quantities H r and at a given axial coordinate z into variables V and I 
localized to a single "point" (the corresponding z coordinate) on a 
Hankel space transmission line. It then takes an infinite number of 
transmission lines to account for the infinite number of possible values 
of the radius in the real space problem. 

To derive an expression for the characteristic impedance of one of 
these Hankel space transmission lines, consider the case of a line having 
only a single frequency carp lex sinusoid traveling in the direction of 
increasing z. In agreement with ordinary transmission line theory, the 
negative sign in the exponents in equations (11) and (12) denotes propa- 
gation in the direction of increasing z. Then equations (11) and (12) 
become 

•A. a - tfa 
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V = e OH) 


Substitute the expressions for the current and voltage from equations 
(13) and (14) into equation (7) 


A 


-I 


o+ 


-Xa 

Ve. 


(<r 



e. 


17 



From this equation, the defining expression for the characteristic imped- 
ance of a Hankel space transmission line is obtained as 



£ 
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This result differs by a negative sign from the expression given in 
El-Markabi and Freeman [4] for the line's characteristic impedance. One 
could argue that the negative branch of the square root function in the 
expression for the propagation coefficient gairma should be chosen, which 
would eliminate the negative sign in the characteristic impedance. How- 
ever, calculations performed with a negative characteristic impedance 
result in predictions that are essentially in agreement with experi- 
mental measurements, whereas the use of a positive characteristic imped- 
ance predicts results that are not in agreement with experiment. Alter- 
natively, one could propose that it is the positive sign that must be 
chosen for the exponents in equations (11) and (12) to correspond to 
propagation in the direction of increasing z. Such an assumption also 
results in predictions that are not in agreement with experiment. 

Transformation of real space current sources into Hankel space is 
discussed in El-Markabi and Freeman [4]. They show that a disc of azimu- 
thal ly directed uniform surface current of phasor value I, with an inner 
radius and an outer radius R 2 » located on the z-axis, becomes a 
sinusoidal current source with phasor value 



connected in parallel with the Hankel space transmission line. The z 
coordinate of this current source is the same as the z coordinate of the 
real space surface current from which it arose. 

Three of these ideal azimuthal ly directed discs of uniform surface 
current have been chosen to model the current distribution in the coil. 
The inner (outer) radius of the disc is equal to the inner (outer) radius 
of the coil. The middle disc is located at the axial center of the coil, 
while the two outer discs are each spaced out a distance of one-third of 
the coil's width from the center of the coil. See Figure 3. Clearly, 
some error occurs here due to the localization of the model's current to 
these discs. 



FIGURE 3 

CURRIHT DISC MODEL 
OF COIL 


Figure 4 shows the Hankel space transmission line configuration 
corresponding to the three current sheet model of the coil next to the 
target. All of the Hankel space calculations are based upon this model. 
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CHAPTER POUR 


ANALYSIS STRUCTURE 


4.1 INTRODUCTION 

In this chapter, the transmission line model of the coil and target 
developed in Chapter 3 and shown in Figure 4 is examined mathematically 
using conventional transmission line theory to predict the Fourier-Hankel 
space behavior of the model when it is in the sinusoidal steady state. 

The analyses performed in this chapter will assume each of the three cur- 
rent sources in Figure 4 has unit phasor value. In view of the analogy 
developed in Chapter 3, use will be made of the terms current and voltage 
in connection with the transmission line model to stand for the Fourier- 
Hankel transforms of the radial magnetic induction Bj-Cz.r.t) and the 
azimuthal electric intensity E^z.r.t) respectively. Throughout this 
chapter, reference should be made to Figure 4 for insight into the equa- 
tions written. 

4.2 METHOD OF CALCULATING COIL IMPEDANCE 

Although the part of the coil impedance that is independent of the 
properties of the ribbon conductor used to wind the coil is calculated in 
Hankel space, as shown in this section, one must begin the derivation of 
the expresion for the impedance in Fourier space. Only in the frequency 
domain does the concept of impedance make sense. The coil impedance will 
be calculated at a frequency omega by exciting the coil with a current 
source of phasor value 1, and finding the phasor voltage response at the 
coil terminals. The coil impedance is then numerically equal to this 
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voltage response. 

The real (Fourier) space model of the coil is three ideal discs of 
azimuthal ly directed uniform surface current. The electric fields in- 
duced in the different parts of this model need to be related to the 
voltage that exists between the coil terminals. To do this, the voltage 
appearing at the coil terminals will be approximated by the arithmetic 
average of the voltages induced at the locations of the three discs of 
the model. Symbolically, 


A/ 






, V ^ /V/ 

i ( - V, 4 . V A ) 


f \ / 

where represents the voltage induced at the location of the disc with 
axial coordinate z^. 

Because of the symmetry of the model, the voltage induced on a cir- 
cular path (centered on the z axis) of radius r and axial coordinate z k 
will be given by 


2 vr (z k , r, <xi) 

The average voltage induced on the infinite number of circular paths 
between r=R^ and r=R 2 i- s 

rT~Ht I H ^' rrr E 4> (2 k 1 r ; <°) 

A » R. ( 

Multiplication of this expression by N, the number of turns in the actual 
coil, yields 


A/ 

V,. 


SL'rt NJ 

R,- 


$ * r E^(z kl r, e*>) dr 


which is the model-predicted voltage induced on an infinitesimally thin N 
turn coil having an inner radius and an outer radius R 2 located at z k . 
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Averaging these induced voltages over the three discs results in the 
phasor coil voltage 


V 



■2-rrKi 

R.-R, 



A 

l- o 


(». , r. oj 


0 Jr} 


V = 


£ 


(W) 


Because this expression is numerically equal to the coil impedance, the 
symbol V will be replaced by the symbol Z in further appearances of this 
expression. The E fields that appear in (15) are frequency domain (Four- 
ier space) E fields. In terms of the Fourier-Hankel space E fields, 
these Fourier space E fields can be written as 


r°° i 

> r > 40 ) 3 \ E.4, (a^ , X , ui) X T, (Xr) dU (ifc) 

O 

where again the list of arguments is used to distinguish between differ- 
ent functions. Specifically, E^z^r,^ ) is the Fourier transform of the 
real space azimuthal electric intensity E^z^^t), and E^z^, X,<o) is 
the Fourier-Hankel transform of the real space azimuthal electric inten- 
sity E$(zk» r » t )* Substituting (16) into (15) yields 


Z 3 


3 


Air hi 


| J E- 4 (b. , X,a>) X 3" ( (Xr) dx] rJr 


Interchanging the order of integration. 


wW ^ ^ 

2 - " [ 3 (R -R j f r dr] £ X dX 

w o < u R i-o 

00 3 

2 - a-K J K (x) £ (*i , X , to) x Jx ( 17 ) 


where 
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K ' U) * \* T i M n ^ 

is the Hankel transformed current corresponding to an infinitesimally 
thin N turn coil of inner radius and outer radius R 2 * carrying a 
phasor current of strength 1/3. 

Expression (17) contains three Fourier-Hankel space electric fields 
(transmission line voltages) that conceptually arose from Fourier space 
phasor currents of value 1/3 on each of the three discs in Figure 3. The 
"scaling factor" that relates Fourier-Hankel space fields due to unity 
Fourier space phasor coil current to Fourier-Hankel space voltages and 
currents due to unity phasor current in each of the three current sources 
in Figure 4 is K'( A ). Then in terms of the Fourier-Hankel space voltage 
^tek* X , co ) due to unity phasor current in each of the Figure 4 current 
sources, (17) becomes 

Z - a-rr f [ K'(x)] £ Z'f (*i, X, co) X c*A 08) 

0 O 1.30 

Expression (18) is the result that will be used to calculate the terminal 
impedance of the coil. 

4.3 RADIAL MAGNETIC INDUCTION AND AZIMUTHAL ELECTRIC INTENSITY 
4.3.1 Target Surface Facing Coil 

Using the results derived in Chapter 3, the characteristic impedance 
and the complex propagation constant r»f the air transmission line and the 
metal transmission line, referred to hereafter : the air line and the 
metal line respectively, will be calculated. For the air line, the 
conductivity sigma is equal to zero, and so 
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Z 


a 


- - 

V^+ 


“ -/X 3 +■ J* 0 /* 0- “ * 

For the metal line, no simplification is possible and so 

-z _ , 

/t m - -fx* + j u y A<r 

^ r*i ~ 4- £«yMr 

According to the analogy developed in Chapter 3, calculation of the 
total current and voltage at z=0 in Figure 4 is equivalent to calculation 
of the radial magnetic induction Bj- and the azimuthal electric intensity 
Eq on the coil side face of the target. Ihis current and voltage calcu- 
lation is performed by replacing the transmission line configuration for 
all z>0 with its equivalent input impedance. Note that both sides of the 
metal line are connected to infinite lengths of air transmission line 
with characteristic impedance Z a . The equivalent impedance looking to 
the right at z=0 in Figure 4 is 


Z <o) 




Zg + Z m ton la (y,v,d) 


(19) 


Using this equivalent impedance, the current reflection coefficient 
at z=0 as seen from the air line is 


<?h 


2L« - 


H a + Z(o) 

and the total current at z=0 is then 


iwi - line + <*» 

in terms of the incident current ti nc , which must be due solely to the 
three current sources. Expression (20) may be simplified as follows. 
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H a - 1 (o) V f 

~Z a + 2(o) ' ■*•*«« 


3 Zg 
+ 2(0) 


A 
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Substituting (19) into (21) yields 


•A, 

^+o+ql 


5 


A 


I. 




Za 

2 + Z S-a* (¥ w cQ 

" Z m + Z a W, (*„<*) 


which, upon replacing Z a and 2^ by their defining expressions gives 


A 

+0+al 


31, 


IrtC 


I + 


Vxh 


_ 1 + jgjjgg t^hQi^J) 


(«) 


Now consider the current source in Figure 4 closest to the target. 
This source produces an incident current on the metal line at z=0 given 
by 


O 



where the factor 1/2 canes from the equal division of the unity phasor 
amplitude current to both sides of the air line containing the current 
source. Tof?i cnr"ent incident at z=0 is by superposition given by 
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(33) 


Substituting (23) into (22) yields 
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for the total steady state current Itotal at 2=0 due to unit strength 
coup lex sinusoidal excitation in each of the three current sources. Mul- 
tiplying (24) by k'(A)I(co), where I(oj) is the Fourier transform of the 
current in the coil (and multiplication by K'( X ) transforms a unity 
Fourier space current into the corresponding Fourier-Hankel space cur- 
rent), yields the desired result 


I(z*o) 


-\q . ~hh 

IM K'M e (» + e 3 + 


x 

I + x 1 * Vx a + joyAO- 
Yx , + j <u / Aflr A + 

Yx j + ioyua- 


- aXb 

e- 3 ) 

+anh (tf m el) 


+<w>li (* m d) 


(a?) 


for the Fourier-Hankel transform of the radial magnetic induction Bj- on 
the coil side face of the metal target when the coil current spectrum is 
I(o>). 

Calculation of the voltage is almost identical to the calculation of 
the current performed above. The voltage reflection coefficient is 
the negative of the current reflection coefficient, 

£(o) - 

~ Z(o) + z a • 

Total voltage at z=0 is 
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in terms of the incident voltage Vi nc . Expression (26) may be simplified 
as 


N/ +o+«| 



2(o)- l a \ Q 

Z(o) + Z a / "* c 


az(o) 

2 a + 2(o) 
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Substituting (19) into (27), 
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V 
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2g + fault (V m «l) 
2m Z a +artta ( d) 


^ + z. 
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+ 2 a +«nVi ( *„,<*) 


which, upon replacing Z a and by their defining expressions, gives 
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Total incident voltage V^ nc at z=0 is given by 
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using (23) for Ij nc . Substituting this expression for Vj nc into (28) 
yields 

\ - Xk _ alV> 
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for the total voltage at z=0 due to unit strength complex sinusoidal 
excitation in each of the three current sources. Multiplying (29) by 
K'(X )!(<*>) yields the desired result 


V ( « * o) 


\<oa \q -^1 - 

XMKVM £• 3 0 + e. 3 •* e. 3 ) 

l + -/x a + + tanh (^d) 

X I + 


( 30 ) 


for the Fourier-Hankel transform of the azimuthal electric intensity ] 
on the coil side face of the target when the coil current spectrum is 
I( o > ). 


4.3.2 Target Surface Opposite Coil 

Knowing the total voltage and current at z=0 (derived in Section 
4.3.1 and given by (30) and (25) respectively) allows a simple trans- 
mission line inverse chain matrix calculation of the total voltage and 
current at z=d, which are the Fourier-Hankel transforms of the azimuthal 
electric intensity and the radial magnetic induction Bj. on the surface 
of the metal target opposite the coil. Denoting these transforms of the 

A 

fields on the coil side face of the target by V n (given by (30) above) 

A 

and I n (given by (25) above). 


V A cosh ( 

- X„ si-h (V„J) 

(31) 

() sinh 

+ i„ <M) 

(35) 
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where Vf and If denote total voltage (electric intensity) and current 
(magnetic induction) on the surface of the target opposite the coil. 

4.4 CALCULATION OF THE AXIAL MAGNETIC FIELD USING THE 

TRANSMISSION LINE MODEL 

The procedure used for calculating the force on the target is de- 
scribed in the next section, and requires both the radial and the axial 
components of the magnetic induction on the target surface next to the 
coil. The radial component of the magnetic intensity is available from 
the transmission line model by performing inverse Hankel and Fourier 
transformations on the calculated transmission line current at z=0, as 
discussed in Section 4.3.1. By using equation (2) in Chapter 3, it is 
also possible to calculate the axial component of the magnetic intensity 
(which does not have a transmission line analog) from the transmission 
line voltage. This will now be shown. 

The azimuthal electric field is given in Fourier space by 

E^r,^) = tf'{V (»,*,*>)} X £v(*>,co) X J t (Xr) olx (33) 

Substituting (33) into equation (2) , 

Jr r b { r 

■pp h { r x W ^ } nu 

Interchanging the order of partial differentiation and integration in 
(34), 

= jsp T,(Xr)} V(,,).,co) (35) 
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Using the identity 


^ ~ x J o^ x ) 

(35) becomes 

W z ( z i r ,^) = ~JJ j V(*,X,a>) \ A T 0 (Xr) dX (3&) 

Equation (36) states that the axial component of the magnetic intensity 
H z in Fourier space is given by the zero order inverse Hankel transform 
of the transmission line voltage at the corresponding z coordinate multi- 
plied by lambda, with the result divided by . Inverse Fourier 
transformation of the right hand side of (36) yields the desired time 
domain axial magnetic intensity. 


4.5 FORCE BETWEEN TARGET AND COIL 

The procedure for calculating the total force between the target and 
the coil utilizes the Maxwell stress tensor, and is performed in real 
space (using the time domain fields). Stratton [41] shows that the total 
force F transmitted by a time varying electromagnetic field across a 
closed surface S is given by 

F = £ [ £(£«>£)£ + ^*(8 o rf) B - + ] da 

where n is a unit vector normal to the surface. The closed surface is 
taken to be the plane z=0 (the coil side face of the target, "closed" at 
infinity). Since n=z, this integral reduces to 

oo Alt 

p = I f 
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The total force tending to separate coil and target is just the z compo- 
nent of this force. 




r ( !*. - ai 

J v au a 

r ao / 


i^£- j r dr 


F* ~ jr J ( - B* ) r dr “ ^ E.J r- Jr ( 37 ) 

Calculations of the fields E^, Bj-, and B z in the prototype experi- 
mental EIDI configuration described in Chapter 5 showed that the second 
integral in (37) is insignificant compared to the first integral. Some 
feeling for why this is true can be obtained by comparing the constants 
that multiply each of the integrals in (37). 
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Accordingly, (37) may be approximated by 

F m <*) * - 5 * J [&,*(*) - B*(*)] r dr ( 38 ) 

Equation (38) is the result used to calculate total force versus time. 


4.6 IMPULS E DEL IVERED 10 TARGET 

The irtpulse delivered to the stationary target is by definition 
given by the integral 

T = At (39) 

J O 

where F z (t) is calculated using (38). 
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CHAPTER FIVE 


A SPECIFIC SYSTEM EXAMPLE, INCLUDING EXPERIMENTAL RESULTS 
5.1 DEFINITION OF THE SYSTEM 

The prototype EIDI system constructed at The Wichita State Uni- 
versity, and the results of the tests made on that system, have been 
described in detail by Dr. Robert Schrag [3]. Most of the material in 
this Chapter has been excerpted from Dr. Schrag's paper. 

Figure 5 shows the prototype EIDI energy discharge system, omitting 
the capacitor charging circuit and the thyristor firing circuit. Two 
identical pulsing coils were operated in series, because that was the 
arrangement used in most of the de-icing tests. However, only one of the 
two coils was utilized in the coil-target assembly, which is detailed in 
Figure 6. 

The effective gap between the coil (copper) surface and the near 
surface of the target was .078". A .032" thick 2024 T3 Aluminum disc was 
used as the target. Diameter of this disc was 5 inches. Two .05" thick 
phenolic spacer plates were used, one to maintain a fixed distance be- 
tween the coil and the target, and the other to maintain the distance 
between the target and the rigid wooden support that prevented motion of 
the target. These plates could be removed, and a special magnetic field 
measuring plate (described in Section 5.4) inserted in their place in 
order to make measurements of the magnetic induction close to the surface 
of the target. Each coil consisted of 30 turns of .024" X .188" rectan- 
gular copper wire spirally wound in a single layer from an inner radius 
of .125" to an outer radius of 1". The initial capacitor voltage 
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utilized for the experimental study was 400 volts. 


5.2 COIL IMPEDANCE MEASUREMENTS 

Impedance measurements on the coil next to the target were made 
using an inpedance bridge. Hie inductance determined by these mea- 
surements, for several frequencies, appears in the table below. Inped- 
ance measurements were also made on the coil when the metal target was 
removed. The real part of the inpedance measured on the coil without the 
target in place was subtracted from the real part of the inpedance mea- 
sured on the coil with the target in place. This increase in resistance 
due to the target is given in the table below. 

RESISTANCE INCREASE AND INDUCTANCE - COIL AND METAL TARGET 

Frequency Inductance Resistance 

(Hertz) (Microhenries) Increase (Milliohms) 


500 

18.6 

8.2 

1000 

17.0 

23.0 

2000 

14.6 

48.0 

4000 

12.6 

77.0 


5.3 CURRENT WAVEFORM 

Initially, current in the coil was measured indirectly by measuring 
the voltage across the .001 ohm non-inductive resistor in Figure 5. 
Difficulties with this approach prompted the purchase of a current trans- 
former from Pearson Electronics, Inc. Using this transformer and a stor- 
age oscilloscope, the current shown in Figure 7 was observed. 
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FIGURE 7 
COIL CURRENT 


5.4 MAGNETIC FTEID MEASUREMENTS 

A magnetic field measuring plate was constructed in the manner 
illustrated in Figure 8. Shallow concentric grooves were cut into both 
sides of a .05" phenolic disc, with radius increments of .2", starting 


TYPICAL GROOVES 
WITH SINGLE 

turn wire 

LOOPS 


CHANNELS WITH 
TWISTED LEADS 



7 DIAMETER 
PHENOLIC DISC, 
.OS" THICK 


SOLDER TAR.S 


FIGURE 8 

PARTIAL ILLUSTRATION 
FIELD MEASURING PLATE 


at r=.2" and ending at r=2.0". Single turn loops of .006" diameter wire 
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were then cemented into these grooves, and their twisted leads brought 
out to solder tabs through radial channels. 

For measuring the fields on either side of the target, the measuring 
plate simply substituted for the corresponding phenolic spacer plate in 
Figure 6. A measurement of the axial flux density was derived from the 
induced voltage in any two neighboring loops connected in series opposi- 
tion. For the two loops illustrated in Figure 8, for example, 

ir('t) d'H 

a /rr(r - r- ) 

v a. i ' 

where B z is in teslas, r is in inches, and u* is in volts. This value is 
the average axial flux density over the area between the two induction 
loops. In the further use of this result, this flux density will be 
assumed to apply at a radius midway between the two loops. 

To measure the radial component of flux density at any radius, the 
front and back loops at that radius are connected in series opposition, 
and calculations are made from 
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where r is the radius of the two induction loops and h is their separa- 
tion, both in inches. 

Plots of the magnetic induction fields obtained from these tests 
with the target in place are shown at selected radii in Figures 9, 10, 

11, and 12. All Bj- data showed an anomolous behavior (irregularities) 
at r=.4" relative to r=.6". A separate check was made, in which the 
plate was reversed (interchanging the two sides). This produced a rever- 
sal of the irregularities, so the effect was probably due to an 
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FAR SIDE OF TARGET FAR SIDE OF TARGET 



inaccuracy in the construction of the plate. 

5.5 MEASUREMENT OF IMPULSE TO TARGET 

The final step in the experiment was the measurement of the inpulse 
delivered to the target when the capacitor was discharged through the 
coil. This was done indirectly, using a ballistic pendulum. The inpulse 
so measured was approximately .012 lb-sec. 
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CHAPTER SIX 


COMPUTER ANALYSIS ON THE SYSTEM EXAMPLE 

6.1 INTRODUCTION 

Theoretical background for the numerical computations described in 
this chapter were developed in Chapters 3 and 4. This chapter will con- 
centrate on a description of the numerical methods used to implement the 
theoretical development contained in those earlier chapters to predict 
the performance of the prototype EIDI system described in Chapter 5. The 
sections that follow are arranged in the order of the Figure 1 Analysis 
Flow Diagram, beginning with Block 4 of that diagram. This is the order 
in which the computations were actually accomplished. 

All computations were performed in FORTRAN 66 on an IBM 370. Source 
code was compiled with the IBM furnished G level compiler. 

6.2 COIL IMPEDANCE 

Equation (18) of Chapter 4, repeated below as equation (40), is the 
coil impedance predicted by the transmission line model. Numerical 
evaluation of the integral in (40) is discussed in this section. 

oo 

ZN = 2 * j" [K'(x)] < I ^ M 

o c * o 

Two common methods are in use to accomplish quadrature of an im- 
proper integral such as the integral above [21]. In the first method, a 
transformation of variables is performed prior to construction of an 
algorithm for estimating the value of the integral. This transformation 
is chosen in such a manner that the new integral has finite limits. The 
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second method simply replaces the infinite upper limit with a finite 
upper limit, selected such that the part of the integral thus ignored 
contributes little to the true value of the integral. This method can be 
employed only if the integrand decays sufficiently rapidly as the vari- 
able of integration increases. Since the controlling factor in the 
leading behavior of the asymptotic expansion, as lambda tends to infini- 

-\q 

ty, of the integrand in (40) is e. , with g a constant, the second 
method was chosen for use in the numerical evaluation of (40). 

A large amount of high quality mathematical software is available 
today [22], [23]. Because of the complexity and cost of writing quality 
mathematical algorithms, most numerical analysts suggest that complex 
scientific calculations be performed using algorithms written by experts 
[24], [25]. An 8 panel adaptive Newton-Cotes algorithm entitled QUANC8, 
described in [25], was chosen to perform the quadrature in (40). A user 
of QUANC8 may select the relative and absolute error performances de- 
sired, and the program then attempts to estimate the value of the inte- 
gral within the selected error criteria. One of the parameters in the 
subroutine QUANC8 is an output variable that contains an estimated error 
bound on the returned value. 

Direct computer evaluation of (40) consumes a large amount of CPU 
time, and is consequently expensive. This is due mostly to the ap- 
pearance of the factor [K'( X )] 2 in the integrand. (Note that K'(X) is 
defined by an integral containing a Bessel function. This makes K' osc- 
illatory, shown in Figure 13.) It was possible to decrease this cost 
considerably by calculating K'( X ) using QUANC8 and fitting a cubic 
spline function to the calculated k'(A) for use in evaluating (40). The 
algorithms used for generating the spline function coefficients and for 
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evaluating the spline function at a given argument are entitled SPLINE 
and SEVAL, respectively, and are described in [25]. 



Once a spline function approximation for k'(X) was available, 

QUANC8 was used to individually estimate the real part (the coil resist- 
ance increase due to the metal target) and the imaginary part (the coil 
reactance) of (40) at selected frequencies. Initial estimates of the 
real and imaginary parts of the integral in (40) were obtained using an 
upper limit of 680. Following this, the programs were run again with an 
upper limit of 240. Little difference was seen in the results of the two 
calculations, leading to the conclusion that the coil resistance and 
inductance estimates were good. Figure 14 shows the effects of the upper 
limit of integration in evaluating the imaginary part of the integral in 
(40) for three different frequencies. 

Figures 15 and 16 show the resistance increase and inductance calcu- 
lated as described above. Measured values of these parameters are also 
shown for comparison. 
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FIGURE 14 

CALCULATED INDUCTANCE 
VS. 

UPPER INTEGRATION LIMIT 


6.3 CURRENT WAVEFORM 

6.3.1 Introduction 

As discussed in Chapter 2, three factors preclude a simple calcu- 
lation of the coil current in the EIDI prototype system. These three 
factors are the nonlinear diode and SCR, and the presence of the metal 
target next to the coil. 

The effects of the target on the coil were taken into account by 
calculating the coil resistance increase and inductance at several fre- 
quencies, as described in the previous section. These two frequency de- 
pendent parameters were then approximated with cubic spline functions, 
again using the subroutine SPLINE. This provided the ability to calcu- 
late computational ly inexpensive yet sufficiently accurate coil imped- 
ances for use in frequency domain circuit calculations. 

6.3.2 Current Before Clamp Diode Conduction 

Although the circuit is nonlinear due to the diode and SCR, it is 
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FIGURE 15 

OOIL RESISTANCE INCREASE DUE TO METAL TARGET 



FIGURE 16 

OOIL INDUCTANCE IN PRESENCE OF METAL TARGET 
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amenable to treatment in a piecewise linear fashion. The frequency 
domain model of the circuit, valid between the time the SCR is initially 
triggered and the time the capacitor voltage becomes zero, is shown in 
Figure 17. In this Figure, the EIDI system capacitor is modeled as an 
ideal discharged capacitor in series with a voltage source. This voltage 





I Uo>) 

FIGURE 17 
FREQUENCY DOMAIN 
CIRCUIT MODEL 

ju>Uui) CLAMP DIODE OFF 


source has a value of zero volts for all time prior to t=0. At t=0, it 
instantly rises to the amplitude and polarity V Q of the initial voltage 
on the actual capacitor, modeling the triggering of the SCR into instan- 
taneous full forward conduction. Voltages and currents calculated from 
this model are good approximations to their corresponding quantities in 
the physical EIDI prototype circuit as long as the diode in parallel with 
the capacitor is not conducting. 

The frequency dependent resistor appearing in Figure 17 is a com- 
posite lumped model of several energy loss mechanisms in the physical 
circuit. It includes the loss in the coil due to the presence of the 
metal target next to the coil (calculated as described in the previous 
section), resistance of the ribbon conductor used to wind the coil (cor- 
rected for skin effect), and the resistance of the cable used to connect 
the coil and capacitor. Resistance in the cable was modeled as frequency 
independent, with a measured D.C. value of .054 ohms. Coil winding 
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resistance was calculated from 


r ac = r dc if* -*(*4*;} 

where 

r DC = D * c * c °il resistance = .0235 ohms 
h = coil thickness = .00477 meters 
S = copper skin depth = I /-/ntTar 
f = frequency in Hertz 

= copper permeability = 4-rf X 10 -7 henries/meter 
<r = copper conductivity = 3.48 X 10 7 mhos/meter 

Resistance R AC is present in both the coil next to the target and the 

idler coil (see Figure 5). 

Similarly, the frequency dependent inductance L in Figure 17 arises 
from more than one source. First, there is the inductance of the coil 
next to the target (calculated in the section above). Second, the idler 
coil and cable connecting the capacitor to the coil had a combined mea- 
sured inductance of 23 microhenries, which was assumed to be independent 
of frequency. As these two inductances are in series in the circuit with 
negligible mutual inductance, they are added to obtain a value for L. 

The leakage reactance and measuring circuit impedance reflected into 
the actual circuit by the current transformer used to measure the current 
were felt to be too small to be significant, and were not included in the 
Figure 17 model. 

By inspection, the Fourier transform of the current in Figure 17 is 
given by 

Vo £(<*>) + ] 

pc + 

cv 0 L 2 L IlL ^! dL * *3 

l + j<oCR(<o) - a^CL(^) 


IM 

IM 
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Note that the delta distribution multiplies its own argument, and conse- 
quently contributes nothing to the time domain current i(t). Equation 
(41) is the basis for the current calculation, valid until the capacitor 
voltage becomes zero and the clamp diode in parallel with the capacitor 
begins conducting. For simplicity and cost considerations, an inverse 
discrete Fourier transform (IDFT) [36] was chosen to approximately invert 
I(o>). To sufficiently minimize the effects of the IDFT approximation to 
the continuous inverse Fourier transform, a 1024 point transform with a 
sample time of 5 microseconds was chosen. This makes the folding fre- 
quency 100kHz, which is considerably above any significant frequencies 
measured in the spectrum of the current in the prototype EIDI system. 

The 200kHz frequency domain window of the IDFT is sufficiently large that 
"windowing" effect errors in the time domain response are not too great. 
(These errors are primarily manifested as Gibbs' ripples on the intitial 
current calculated in Section 6.3.3, described in Section 6.3.4 and shown 
in Figure 19.) 

In order to determine how long a time this calculated current ap- 
proximates the actual circuit current, a second calculation was performed 
using the circuit model in Figure 17. The model predicts a frequency do- 
main voltage corresponding to the physical EIDI capacitor voltage of 
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The term involving the delta distribution again contributes nothing to 
the inverse Fourier integral yielding the time domain voltage, and is 
dropped from consideration to yield (42). A 1024 point 5 microsecond 
sample time IDFT was used to approximately calculate the inverse Fourier 
transform of (42). Output from the IDFT showed that the capacitor volt- 
age would slew negatively through zero volts at t=279 microseconds. This 
is the time at which the diode across the capacitor is modeled as coming 
into conduction and acting as a short circuit. Beyond this time the 
model of Figure 17 is no longer valid, and consequently the current 
predicted by the model is no longer correct. 


6.3.3 Cu r rent After Clanp Diode Conduction 

Once the diode comes into conduction, the model of the EIDI circuit 
for all future time is as shown in Figure 18. The inductance and resis- 





FIGURE 18 
FREQUENCY DOMAIN 
CIRCUIT MODEL 
CLAMP DIODE ON 


tance in this circuit have the same physical significance that they had 
in Figure 17, and their values are calculated for any desired frequency 
using the same algorithms. For computational convenience, time is "re- 
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set" to zero in this circuit, even though the circuit does not come into 
existence until t=279 microseconds in the Figure 17 circuit. Since the 
current in the coil is not initially zero, the circuit's magnetic energy 
storage at the new "t=0" is modeled by including a DC current source in 
parallel with a frequency dependent inductance. This current source is 
zero for all time prior to t=0, and for all future time has a D.C. value 
I Q equal to the value of the current in the Figure 17 circuit at the time 
the capacitor voltage became zero. 

Current in the Figure 18 circuit is given in the frequency domain by 
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Note that the term containing the delta distribution once again con- 
tributes nothing to the inverse Fourier transform of I(oo ), and is drop- 
ped in (43). An IDFT was used to calculate the approximate inverse Four- 
ier transform of I (co ) given in (43), in exactly the same manner that the 
expression for I(oo) in (42) was inverted. The program is very similar 
to the one used in the Section 6.3.2 current calculation. 


6.3.4 Combining Pre- and Post Clamp Diode Conduction Current 

Because of the presence of the current source in parallel with the 
inductor in Figure 18, the current in this circuit model is discontinuous 
at t=0. Consequently, the current cannot be bandlimited. An inherent 
assumption in the use of the discrete Fourier transform is that the 
signal is bandlimited. The time domain current calculated by the IDFT 
showed a small amount of ripple during the first 100 microseconds due to 
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the current spectrum not being bandlimited. Prior to "joining" in time 
the current predicted by the Figure 17 model with the current predicted 
by the Figure 18 model, this ripple was eliminated by graphically choos- 
ing current values lying close to the IDFT predicted values that joined 
smoothly with the current from the Figure 17 model. Figure 19 illus- 
trates this procedure. Current ripple amplitude was 54 amps peak to peak 
at t=10 microseconds, decaying to 1 amp peak to peak at t=75 microsec- 
onds, with time measured in the Figure 18 model. Figure 20 shows how 
good the agreement is between this piecewise linear model predicted 
current and the current measured in the prototype EIDI circuit. 

With the model predicted current now available at 5 microsecond 
intervals from t=0 to t=5.12 milliseconds, a 1024 point 5 microsecond 
sample time DFT was performed on the current samples to estimate the 
Fourier transform I(oa) of the model current. This transform is needed 
for use in calculating the magnetic induction fields from the Hankel 
space transmission line coil-target model, as discussed in the next 
section . 
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PRE- AND POST- CLAMP DIODE CONDUCTION CURRENT 
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FIGURE 20 
COIL CURRENT 


6.4 MAGNETIC INDUCTION 

6.4.1 Introduction 

Knowledge of the radial and axial magnetic induction on the coil 
side face of the target is required for calculation of the total force on 
the target. Although conceptually simple, the numerical calculation of 
these two fields presented more computational difficulties than were 
encountered in the combined total of all other calculations performed. 
Initial attempts at calculating the magnetic induction, performed using 
single precision arithmetic and using professionally written quadrature 
routines, produced results that were incorrect by orders of magnitude. 

6.4.2 Radial Magnetic Induction on the Coil Side Face of the Target 
Expression (44) shows the iterated improper integral that is to be 

evaluated numerically to calculate the time domain radial magnetic induc- 
tion Bj-^r.t). This integral is the inverse Fourier-Hankel transform of 
the total current at z=0 given by (25) of Chapter 4. 
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Quadrature of iterated integrals is almost always difficult [21]. 

One of the most common methods of evaluating such integrals, the Monte 
Carlo method, could not even be considered for use in (44) due to the 
enormity of the very expensive and random complex function evaluations 
required by such an approach. Furthermore, the integrand in (44) is 
extremely oscillatory, having both positive and negative complex values, 
due to the Bessel function in the inverse Hankel transform, the function 
K'( A ), and the complex exponential in the inverse Fourier transform. 
Quadrature of such integrands is all but impossible using Monte Carlo 
methods due to extreme smearing [25], [26], It is the oscillatory nature 
of the integrand in (44) that makes its evaluation difficult. 

After much trial and error, a workable approach to quadrature of the 
integral in (44) was obtained, and is described in this section. No 
claim is made for optimality or near optimality in this approach. After 
many computations of the magnetic induction had been performed using this 
procedure, it became apparent that simplifications could be made, while 
still obtaining sufficient accuracy. However, these simplifications have 
not been tested, and the original approach to quadrature of the integral 
in (44) will be described. 
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Measurements of the radial magnetic induction near the coil side 
face of the target in the EIDI prototype system provide some insight into 
a suitable numerical procedure for evaluating the integral in (44). The 
Fourier transforms of these observed fields are smooth (higher order 
derivatives with respect to frequency are small) with most of the energy 
confined to relatively low frequencies. This suggested the use of an 
IDFT to numerically perform the inverse Fourier transformation in (44). 
Further impetus for use of the IDFT is given by noting that the DFT 
procedure for calculating the Fourier transform I(<o) of the current, 
described in the preceding section, yielded transform values at equally 
spaced frequencies. These are the appropriate frequencies for calculat- 
ing a 1024 point IDFT, with a 5 microsecond sample time, of the Fourier 
space radial magnetic induction. However, numerical evaluation of the 
inverse Hankel transform in (44) at the 513 frequencies needed for calcu- 
lating a 1024 point IDFT yielding a real sequence is prohibitively expen- 
sive. This expense was avoided by using the previously noted fact that 
the expected Fourier transform of the magnetic induction is smooth, with 
most of the energy confined to low frequencies. Accordingly, the inverse 
Hankel transform in (44) was evaluated at only 77 frequencies, chosen to 
provide a reasonable representation of the behavior of the Fourier trans- 
form of the known radial magnetic induction. Spline function approxima- 
tions were then fitted to the real and imaginary parts of the calculated 
inverse Hankel transforms. The IDFT routine to calculate the desired 
time domain magnetic induction uses these spline functions to form 
Bj-tO.r, to ) at the 1024 frequencies needed. 

With a procedure for evaluating the inverse Fourier integral in (44) 
available, a search for a suitable method for quadrature of the inverse 
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Hankel integral was initiated. The infinite upper limit of this integral 

was simply replaced with a suitably large finite upper limit, due to the 

controlling factor in the leading behavior of the asymptotic expansion of 

-X q 

the integrand being e . With the problem of the infinite upper limit 
gone, the remaining difficulties may be divided into two somewhat over- 
lapping classes. 

Selection of a suitable algorithm for performing the quadrature is 
essential if accurate results are to be obtained. Most of the common 
quadrature algorithms are incapable of dealing with oscillatory integrals 
without some help. Even with the best of help, examples of problems 
where these algorithms completely fail abound. This has given rise to 
highly specialized techniques for quadrature involving oscillatory inte- 
grands [27], [28], [29], [30]. The use of one of these specialized 
algorithms was avoided by writing a double precision Gaussian quadrature 
routine [21], [31]. This routine was then used to perform the inverse 
Hankel integration in (44) over the range 0 < X < 10, then over the 
range 10 < X < 20, etc., stopping when the desired upper limit of 
integration had been reached. The results of these integrations, which 
can be interpreted as terms in a sequence that are to be summed, were 
then added together first using straightforward addition, and then using 
the Euler series summation convergence acceleration algorithm DTEUL from 
the NASA Lewis Research Analysis Center Software Library [34]. For each 
radius at which the radial magnetic induction was evaluated, and for each 
of the 77 Fourier frequencies, good agreement was achieved between the 
results of these two summation methods. 

Several different finite upper limits were used to take the place of 
the infinite upper limit in the inverse Hankel quadrature. It was 
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experimentally determined that upper limits greater than 1000 resulted in 
very little change in the calculated induction fields. 

The second problem concerned the precision of the FORTRAN implement- 
ed on the IBM 370. Single precision floating point word length on the 
370 is approximately 7 decimal digits (24 bit mantissa) [25]. This is 
insufficient for nearly all complex scientific calculations [32], [33]. 
Double precision floating point word length on the 370 is approximately 
15 decimal digits (56 bit mantissa) [25], and was used for all computa- 
tions involving the inverse Hankel integral. Without the use of double 
precision arithmetic, inverse Hankel quadrature was impossible due to 
smearing. 

Although the use of double precision greatly reduced the effects of 
smearing, it created problems of its own. A double precision Bessel 
function, and a double precision algorithm for K'(\) become necessary to 
evaluate (44). Simply converting a single precision algorithm for k'(A) 
to double precision does not yield sufficient accuracy to obtain good 
results. Double precision Bessel function algorithms were unavailable at 
The Wichita State University. Double precision algorithms were written 
for J 0 (x) and J^(x). The double precision algorithm that was written for 
K'(X ) is described in the Appendix. Figure 21 provides a comparison 
between the measured radial magnetic induction close to the coil side 
face of the target, and the predicted radial magnetic induction on the 
coil side face of the target. 

6.4.3 Axial Magnetic Induction on the Coil Side Face of the Target 

The magnetic induction B z on the target next to the coil is calcu- 
lated by performing an inverse Fourier transform on (36) in Section 4.4 
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FIGURE 21 
Bj- VERSUS TIME 
NEAR SIDE OF TARGET 


57 



of Chapter 4. Then 
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where V(0, X , u>) is given by (30) in Section 4.3.1 of Chapter 4. Sub- 
stituting (30) into the above. 
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The similarity between (44) and (45) is evident, even though the order of 
the inverse Hankel transforms is different. Because of this similarity, 
the algorithms developed to perform the quadrature in (44) were used, 
with only evident minor changes necessary, to perform the quadrature in 
(45). A comparison at several different radii between the measured axial 
magnetic induction close to the coil side face of the target and the 
calculated axial magnetic induction on the surface of the target closest 
to the coil is shown in Figure 22. 

6.4.4 Magnetic Induction on the Opposite Face of the Target 

Although it was not needed in the calculation of the force between 
the coil and target, the magnetic induction on the side of the target 


58 



© 



FIGURE 22 
B z VERSUS TIME 
NEAR SIDE OF TARGET 
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away from the coil was calculated for comparison with the measured induc- 
tion close to the same target surface. Section 4.3 in Chapter 4 derived 
expressions for the total voltage Vf and total current If at the junction 
between the metal line and air line opposite the current sources. These 
quantities correspond respectively to the azimuthal electric intensity 
and the radial magnetic induction on the surface of the target opposite 
the coil. 

Inverse Fourier-Hankel transformation similar to that in (44), but 
with the integrand corresponding to the current If given in Section 4.3, 
yields the radial magnetic induction. Only minor changes are necessary 
in the FORTRAN program to calculate the induction. Figure 23 provides a 
comparison between the measured and calculated radial magnetic induction 
close to and on the target surface away from the coil. 

Inverse Fourier-Hankel transformation similar to that in (45), but 
with the integrand corresponding to the voltage Vf given in Section 4.3, 
yields the axial magnetic induction. Figure 24 shows the comparison 
between the measured and calculated axial induction close to and on the 
target surface away from the coil. 

6.5 FORCE VERSUS TIME 

It was shown in Section 4.5 of Chapter 4 that the total force be- 
tween the coil and target is given approximately by 

o© 

F * (t) = jr J C&a(o ) r ,t) - B^(o r dr (46) 

For use in this integral, radial magnetic induction was calculated at 5 
microsecond intervals at radii .01", .2", and in .2" increments up to a 
maximum of 2.0", on the surface of the target closest to the coil using 
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FIGURE 23 
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FIGURE 24 
B z VERSUS TIME 
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the algorithms described in Section 6.4.2. Axial magnetic induction was 
also calculated at 5 microsecond intervals at radii .01”, .1”, and in .2” 
increments up to 1.9”, on the same target surface using the algorithms 
described in Section 6.4.3. Cubic spline functions were then fitted to 
the squares of the radial variation of these two magnetic fields at times 
0, 50, 100, 150, .... 900, and 1125, 1400, and 2000 microseconds for use 
in performing the quadrature indicated in (36), using an upper limit of 
1.9". The quadrature was done exactly (to the limit of the machine 
arithmetic) on the spline function approximations to the squared fields 
by integrating the cubic polynomial form of the spline function between 
knots, and using the spline function coefficients in the result. 

Figure 25 shows the force versus time calculated using this proced- 
ure. 


6.6 IMPULSE 

Equation (39) in Chapter 4, repeated below, was used to calculate 

r = C F * (t) dt 

the impulse delivered to the target by the coil. A spline function was 
fitted, with time as the variable, to the force calculated in Section 6.5 
above. An exact integration was performed on the cubic polynomials of 
the spline function between knots, with the result that the inpulse 
calculated was .008 lb-sec. This is lower than the .012 lb-sec inpulse 
measured using a ballistic pendulum. It should be noted that the quadra- 
ture of the force integral in (46) above used an upper limit of r=1.9” 
instead of infinity. Although the induction decays as the radius tends 
to infinity, an unknown part of the forc^e has been ignored by not taking 
the induction fields at radii greater than 1.9" into account. 


63 




CHAPTER SEVEN 


CONCLUSION AND RBCOMflNQAXIONS FOR FURTHER WORK 

7.1 CONCLUSION 

7.1.1 Summary 

A method of modeling the electrical system aspects of a coil and 
metal target configuration resembling a practical Electro-Impulse De- 
Icing installation, and a simple circuit for providing energy to the 
coil, was presented. The model was developed in sufficient theoretical 
detail to allow the generation of computer algorithms for the current in 
the coil, the magnetic induction on both surfaces of the target, the 
force between the coil and target, and the impulse delivered to the 
target. These algorithms were applied to a specific prototype EIDI test 
system for which the current, magnetic fields near the target surfaces, 
and inpulse had previously been measured. 

Coil impedance was the first quantity calculated using the algo- 
rithms. Agreement between the impedance calculated and the impedance 
measured was seen to be very good for the resistive part, and reasonable 
for the reactive part. Despite the simple model used for the circuit 
providing energy to the coil, excellent agreement was obtained between 
the predicted and measured coil current. 

Measured and predicted magnetic induction fields were not directly 
compared, due to the fields having been measured close to, but not on, 
the target surfaces using spatial averaging methods. The only fields 
calculated were on the target surfaces, for use in calculation of the 
force between the coil and target. Nevertheless, there was seen to be 
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very reasonable agreement between measured and calculated magnetic 
fields. The character of the time variation of these fields changed 
considerably with radial distance from the axis, and the algorithms 
correctly predicted these changes. 

Calculation of the impulse easily provided the greatest disagreement 
between a predicted and measured quantity, with a -33% error in the 
calculated impulse. Impulse was measured using a ballistic pendulum 
containing the metal target, so, for this measurement, the metal target 
was no longer held rigid when the capacitor was discharged through the 
coil. This does not satisfy the model assumption of a stationary system. 
However, the period of the pendulum was sufficiently long that, during 
the time interval when most of the force was developed on the target, 
negligible motion of the pendulum should theoretically have occured. 
Motion of the target during the measurement of the impulse is not felt to 
be a satisfactory explanation for the discrepancy between the measured 
and calculated impulse. It was mentioned in Section 6.6, where the 
calculation of the impulse was described, that the infinite upper limit 
in the integral yielding the theoretical impulse had been replaced with a 
finite upper limit for the purpose of quadrature, thus incurring an 
unknown error. This procedure for quadrature of the impulse integral is 
felt to be the most likely source of error in the calculated impulse. 

7.1.2 Original Contributions 

Reference [4] provided most of the basic methods used in this report 
in modeling the interaction between the coil and target. It was shown 
that an error had occured in the definition given in [4] of the charac- 
teristic impedance of a Hankel space transmission line, and the corrected 
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definition was used in the theoretical development of the model. For 
many transmission line calculations, this error does not become evident 
because impedances occur in ratios (e.g., the voltage and current reflec- 
tion coefficient calculations in Section 4.3.1). 

Reference [4] provided no method of calculating the axial magnetic 
induction B z (z,r,t) from the transmission line model. An integral solu- 
tion for this field, in terms of the transmission line voltage, was de- 
rived in Section 4.4. This solution allowed the use of nearly all of the 
algorithms developed for the calculation of the radial magnetic induction 
Bj-(z,r,t) in calculating the axial magnetic induction B z (z,r,t). 

The most significant contribution of this work was the development 
of FORTRAN algorithms for performing the inverse Fourier-Hankel transfor- 
mations yielding the induction fields, described in Section 6.4. While 
there were no theoretical contributions made during the development of 
these algorithms, several diverse results from the fields of numerical 
analysis and approximation theory had to be applied in concert to create 
a working algorithm. During this development, an algorithm for the 
generation of Struve functions of the first and second orders was devel- 
oped that, according to a computer search of the literature, is the most 
accurate reported. 

7.2 RECOMMDTOATIOKS FDR FUTURE WORK 

If the methods developed in this report are to be applicable as 
design tools for EIDI systems, the algorithms for the calculation of the 
magnetic induction fields should be made more efficient. These algo- 
rithms take nearly half an hour of CPU time on the IBM 370 to calculate 
either the axial or the radial magnetic induction versus time at one 
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specific radius. Both the radial and axial components are required, each 
at eleven different radii, for the force calculation described in Section 
6.5. A minimum of eleven hours of CPU time is too great for the evalua- 
tion of the force-time profile of a proposed EIDI configuration. The 
least desirable feature of the methods presented in this report is the 
inordinate CPU time required for the calculation of the induction fields. 
A sophisticated convergence acceleration routine could perhaps be devised 
specifically for more economic calculation of the induction fields. 

Human intervention is required to proceed along the Analysis Flow 
Diagram of Figure 1 (page 4) in evaluating a specific EIDI system. This 
is because the outputs from the various programs, represented by such 
Figure 1 quantities as the coil impedance Z(<o), the current i(t) and 
I ( 60 ) , are not in the form required as the inputs for the program at the 
next block in Figure 1. A significant amount of design automation could 
be accomplished by writing one long program, using as a skeleton the 
programs developed for this report, that will take as input the geometry 
of the coil-target configuration and the circuit used to provide energy 
to the coil, and provide as output a force-time profile and the total 
impulse delivered to the target. 

Not all possibilities have been exhausted in the search for an 
analytical (or semianalytical) solution to the EIDI design problem. 
Perturbation techniques have been suggested as a possible method for 
analytical Hankel inversion in the calculation of the magnetic induction, 
and should be investigated, as this is the area suffering the greatest 
computational expense. Furthermore, analytical solutions (if sufficient- 
ly simple) are often capable of providing insight into a problem not 
provided by numerical solutions. 
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APPENDIX 


CALCULATION OF K 7 ( X ) 


K,(x > - r ,5 L * T|(X)<) 


(M7) 


The function K '( X ) arises from the transformation of real space 
current discs into Hankel space. Using N=30 turns (the number of turns 
of the coil in the experiment described in Chapter 5) in the expression 
for K'( X) on the top of page 24 yields 

R , 

Straightforward quadrature of (47) using the Newton-Cotes algorithm 
QUANC8 was initially performed for generating values of K / ( X) for use in 
approximating K'( X) with a cubic spline function. This provided suffi- 
cient accuracy for use in numerical calculation of coil impedance as de- 
scribed in Section 6.2. Attempts at calculating the radial magnetic in- 
duction Bj- using this spline function approximation for K'( X ) were a 
total failure. This appendix describes a procedure for calculating 
K'( A ) that is accurate to 13 digits when implemented in double precision 
FORTRAN on an IBM 370. 


The integral in (47) can be evaluated in closed form in terms of 
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so that (47) becomes 
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From reference [18], 

I a d a - - m,<»»3-.(3)3 

where H 0 (x) and H^(x) are Struve functions of orders 0 and 1 respectively 
[18], [19], and Jq(x) and J^(x) are Bessel functions of the first kind of 
orders 0 and 1 respectively. Using this result in (48) yields 

h ° ( 8> ' H.WWlL.x* 

^ 1 O l 

In order to use this result, double precision algorithms for generating 
the Bessel and Struve functions must be available. Double precision 
Bessel functions are readily available, but Struve functions are not. A 
computer search of the literature resulted in a reference to a Naval 
Research Laboratory report that contained FORTRAN source code for gene- 
rating integer order Struve functions with positive arguments [36]. This 
code could not be used because it used subroutines to which access was 
not readily available. 

For several years, mathematicians have been aware of the desirabili- 
ty of using truncated Jacobi series of Chebyshev polynomials for numeri- 
cal approximation of various special functions [37], [38]. Luke [37] 
provides coefficients and c n for Chebyshev series expansions 

H.<*) = if b„T,„(i) , 1*1 * 8 

0 = 0 

U,<*> = f c„T^(f) 1*1 i 8 

*» = o 

where T 2 n (x) is a Chebyshev polynomial of the first kind of order 2n. 

All coefficients having magnitudes greater than 10 -20 are given, allowing 
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generation of Hq(x) and H^(x) with 20 decimal digit accuracy for argu- 
ments whose magnitude is less than 8 [39]. Using the identity 

T, „(x) = T„(a* a -0 

and retaining sufficient terms for 15 digit accuracy, these series 
expansions become 

H.M “ if fc.T.Cafi) 4 -'] , IkI s 8 (TO) 

f\ a O 

U, <*> = Z IHI S 8 (51) 

n-o 7 

Then for arguments A such that < 8, the expression on the right 

hand side of (49) is evaluated using (50) and (51). Chebyshev polyno- 
mials are evaluated in double precision arithmetic using the subroutine 
DCNPS from the NASA Lewis Research Analysis Center Software Library [40]. 
Luke also lists coefficients d n and e n for the series 

«.w ' ~ kr E X (i) , w * • 

- O / 

U,(H) - Y,(x) = £ E e,T to (|.) 1*1 i s 

n s o 7 

for 15 digit accuracy. The functions Yq(x) and Y 1 (x) are Bessel func- 
tions of the second kind of orders 0 and 1 respectively. Rather than use 
these series directly to evaluate H 0 (x) and H 1 (x) for arguments greater 
than 8, some simplification is possible. Writing 

X (><) - k EXT,^.!) + Y o (k) = A„(*) + y o (x) 

H.W = £ f *.T Jn (£j * Y,(x) s \<K) + Y,(K) 

and substituting these expressions for H 0 (x) and H^(x) into the expres- 
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sion Ji(x)Hq(x) - H 1 (x)Jq(x) appearing in (49), 

T,(x)M o (k) - U | (x)T > <k) = + Y 0 (x)] - [A.UHY/x)] J o (k) 

= J,Cx) A 6 (x) - A,(x)T 0 (x) + T/x)Y 0 (x) - Y,(x) T 0 (x) 

= T.W A 0 (x) - A,(x) J c (x) + ^ 

using a well known property of the Wronskian of Bessel functions of the 
first and second kinds [18]. Thus, the expression on the right hand side 
of (49) can be written as 

i hVA' ) * V3 )T “ ( i)' + ^ 

when it is to be evaluated at an argument y = XR n > 8. 

Using these results, K'( X ) can be evaluated from (49) for non- 
negative values of the argument to at least 13 digit accuracy. 
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